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Free-Surface  Flow  and 
Fluid-Object  Interaction 
Modeling  With  Emphasis 
on  Ship  Hydrodynamics 

This  paper  presents  our  approach  for  the  computation  of  free-surfacelrigid-hody  inter¬ 
action  phenomena  with  emphasis  on  ship  hydrodynamics.  We  adopt  the  level  set 
approach  to  capture  the  free-surface.  The  rigid  body  is  described  using  six-degree-of- 
freedom  equations  of  motion.  An  interface-tracking  method  is  used  to  handle  the 
interface  between  the  moving  rigid  body  and  the  fluid  domain.  An  Arbitrary 
Lagrangian-Eulerian  version  of  the  residual-based  variational  multiscale  formulation 
for  the  Navier-Stokes  and  level  set  equations  is  employed  in  order  to  accommodate  the 
fluid  domain  motion.  The  free-surface/ rigid  body  problem  is  formulated  and  solved  in  a 
fully  coupled  fashion.  The  numerical  results  illustrate  the  accuracy  and  robustness  of 
the  proposed  approach.  [DOI:  10.1115/1.4005072] 


1  Introduction 

Accurate  prediction  of  wave  loading  on  vessels  and  their  motion 
motions  necessitates  the  development  of  a  simulation  framework 
that  involves  free-surface  flow  and  fluid-structure  interaction  (FSI) 
(see,  for  example.  Ref.  [1]).  The  simulation  framework  must  be 
such  that  it  is  able  to  keep  track  of  two  types  of  interfaces:  the  air- 
water  interface  and  the  fluid-structure  interface. 

Depending  on  flow  conditions,  the  free-surface  motion  may  be 
smooth  or  violent,  with  wave  breaking  and  other  topological 
changes.  As  a  result,  the  use  of  an  interface-capturing  method  (see 
Ref.  [2]  for  the  terminology)  is  a  convenient  and  practical  choice 
for  the  proposed  application.  In  this  work,  we  make  use  of  the 
level-set  method  to  handle  free-surface  flow  (see  Refs.  [3-5]).  In 
the  proposed  methodology,  the  boundary  between  the  two  fluids  is 
described  implicitly  as  a  zero  level  set  of  a  scalar  function  defined 
in  the  problem  domain.  The  subdomains  corresponding  to  nega¬ 
tive  and  positive  values  of  the  level  set  function  are  those  occu¬ 
pied  by  air  and  water,  respectively.  The  level  set  function  is 
simultaneously  a  signed-distance  function,  meaning  its  magnitude 
at  a  point  in  3D  space  is  the  distance  from  that  point  to  the  air- 
water  interface,  and  its  sign  determines  if  the  point  is  in  the  water 
or  air  domain.  The  signed  distance  property  of  the  scalar  function 
is  not  necessary  in  general;  however,  it  has  several  accuracy  bene¬ 
fits,  as  discussed  in  Ref.  [6], 

Simulating  interaction  of  free-surface  flow  with  moving  and  de¬ 
formable  structures  requires  additional  computational  technology 
that  is  able  to  track  the  interface  between  the  structure  and  the  sur¬ 
rounding  air-water  medium.  Therefore,  the  Mixed  Interface- 
Tracking/Interface-Capturing  Technique  (MITICT)  is  used  to 
track  the  fluid-structure  interface  and  to  capture  the  air-water 
interface.  The  MITICT  [7]  was  introduced  primarily  for  fluid- 
object  interactions  with  multiple  fluids.  The  MITICT  was  success- 
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fully  tested  in  Ref.  [8],  where  the  interface-tracking  technique 
used  was  a  space-time  formulation  [9,10],  and  the  interface- 
capturing  method  was  the  Edge-Tracked  Interface  Locator  Tech¬ 
nique  (ETILT)  [7].  It  was  also  tested  in  Ref.  [11]  by  using  a  mov¬ 
ing  Lagrangian  interface  technique  [12]  for  interface  tracking  and 
the  ETILT  for  interface  capturing.  The  interface-tracking  tech¬ 
nique  used  in  the  MITICT  can  also  be  the  Arbitrary  Lagrangian- 
Eulerian  (ALE)  method  [13],  which  is  employed  in  this  work.  In 
this  paper  we  assume  that  the  structures  are  complex-geometry  six 
degree-of-freedom  (6DOF)  rigid  objects.  The  extension  to  flexible 
structures  will  be  pursued  in  future  work. 

The  paper  is  outlined  as  follows.  In  Sec.  2,  we  present  the  gov¬ 
erning  equations  of  free-surface  flow.  In  Sec.  3,  we  present  our 
discrete  free-surface  formulation,  which  includes  the  Residual- 
based  Variational  Multiscale  formulation  (RBVMS)  in  ALE  form, 
weak  enforcement  of  essential  boundary  conditions,  re-distancing, 
and  restoration  of  the  global  fluid  mass  balance.  The  free-surface 
formulation  is  essentially  that  developed  in  Refs.  [14,15].  In  this 
section,  we  also  present  our  rigid  body  dynamics  and  mesh  motion 
formulations.  In  Sec.  4,  we  present  our  time  integration  method 
and  solution  strategy  for  the  coupled  problem.  We  advocate  that 
the  fluid  and  level  set  solutions  are  handled  in  a  fully-coupled 
fashion,  which  significantly  increases  the  robustness  of  the  pro¬ 
posed  methodology.  We  also  present  a  novel  time  integration 
strategy  for  the  rigid  body  equations,  where  the  rotation  matrix 
becomes  an  additional  problem  unknown  that  is  integrated  in  time 
along  with  the  displacement  and  rotational  degrees  of  freedom.  In 
Sec.  5,  we  present  two  numerical  examples.  The  first  one  is  the 
well-known  MARIN  dam  break  problem  [14-18],  which  is  a  free- 
surface  flow  without  fluid-object  interaction.  The  problem  is  used 
to  illustrate  the  effect  of  the  penalty  parameter  in  the  level  set  re¬ 
distancing  equations  and  the  effect  of  using  a  strongly  coupled 
Navier-Stokes/level  set  convection  formulation.  The  second 
example  consists  of  a  DTMB  5415  Navy  combatant  at  lab  scale  in 
head  waves  of  large  amplitude,  which  illustrates  direct 
applicability  of  the  proposed  methodology  to  ship  hydrodynamics 
simulations. 
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2  Governing  Equations  of  the  Free-Surface  Air-Water 
Modeling 

In  this  section,  we  summarize  the  governing  differential  equa¬ 
tions  of  free-surface  flow  on  a  moving  spatial  domain.  Let 
fi,  €  R d ,d  =  2,3  denote  the  combined  air-water  domain  at  cur¬ 
rent  time  t  and  let  T,  be  its  boundary.  The  domain  fi,  is  decom¬ 
posed  into  the  water  and  air  subdomains,  £1“'  and  fi",  respectively, 
and  T™  denotes  the  boundary  between  the  air  and  water  subdo¬ 
mains.  See  Fig.  1  for  an  illustration. 

In  this  work  we  make  use  of  the  level  set  method  to  capture  the 
air-water  interface  (see,  e.g..  Refs.  [3-5]). 

For  this,  we  introduce  a  scalar  function  < p(x,  t)  and  define: 


fi"  =  {x \4>(x,  t )  <0,  Vx  e  fi,} 

(1) 

fi}1  =  {x \4>(x,  t)  >o,  Vx  e  fi,} 

(2) 

r}“-  =  {x| (j)(x,  t)  =  0,  Vx  G  fi,} 

(3) 

In  our  modeling  framework,  the  air  and  water  domains  will  be  dis¬ 
tinguished  by  assigning  the  corresponding  values  to  the  fluid  den¬ 
sity  and  viscosity.  Namely,  we  assume  that  the  density,  p,  of  the 
two-fluid  air-water  medium  is  given  by: 

P  =  P»H(4>)  +pa{l 

(4) 

where  pw  and  pa  are  the  densities  of  water  and  air, 
and  H(4>)  is  the  Heaviside  function  given  by 

respectively, 

/  0  if  4>  <  0 
=  {  1  if  4>  >  0 

(5) 

Likewise,  the  combined  dynamic  viscosity,  p,  is  given  by: 

P  =  PwH(4>)  +  /ta(l  —  H(4>))  (6) 

where  pw  and  pa  are  the  dynamic  viscosities  of  water  and  air, 
respectively. 

With  this  definition  of  the  fluid  material  parameters,  the 
Navier-Stokes  equations  of  incompressible  flow  in  the  arbitrary 
Lagrangian-Eulerian  (ALE)  description  take  the  form  of: 


(  du 

+(u  —  fi)  ■  Vxu  —  f  J  —  V*  ■  a  =  0 

X  / 

(7) 

VA  ■  u  =  0 

(8) 

where  the  fluid  Cauchy  stress,  o,  is  defined  as: 

cr(u,/t)  =  —pi  +  2/(V}u 

(9) 

u  and  p  are  the  fluid  velocity  and  pressure,  f  is  the  body  force  per 
unit  mass,  fi  is  the  velocity  of  the  fluid  domain,  and  Vs  is  a  sym¬ 
metric  gradient.  In  Eq.  (7),  the  partial  time  derivative  is  taken 


Fig.  1  The  fluid  spatial  domain  decomposed  into  the  water 
and  air  subdomains  denoted  by  £2"  and  fif,  respectively.  The 
air-water  interface  is  denoted  by  l  flv. 


with  respect  to  a  referential  coordinate  x  held  fixed.  In  Eqs.  (7) 
and  (8)  the  space  derivatives  are  taken  with  respect  to  the  current 
configuration  spatial  coordinates  denoted  by  x.  We  note  that  the 
fluid  domain  velocity  is  assumed  to  be  completely  independent  of 
the  velocity  of  the  fluid  material  particles. 

We  assume  that  the  air-water  interface  is  moving  with  the 
fluid  material  particles,  which  we  model  by  means  of  an  addi¬ 
tional  convection  equation  for  the  level  set  <j>  in  the  ALE 
description: 


fi)  ■  VX  =  0 


(10) 


The  above  equations,  with  the  associated  boundary  conditions 
constitute  an  air-water  free  surface  formulation  on  a  moving  do¬ 
main  fi,  at  the  continuous  level. 


3  Discrete  Formulation 

3.1  Discrete  Formulation  of  the  Two-Fluid  Problem.  We 

discretize  the  Navier-Stokes  and  level  set  equations  using  the  re¬ 
sidual-based  variational  multi-scale  (RBVMS)  formulation.  The 
RBVMS  formulation  was  originally  proposed  in  Ref.  [19]  and 
applied  to  a  variety  of  situations  involving  fluid  flow  and  fluid- 
structure  interaction  [20-28].  For  a  thorough  derivation  of  the 
RBVMS  formulation  the  reader  is  referred  to  these  publications. 
The  RBVMS  formulation  for  free-surface  flow  was  also  devel¬ 
oped  in  [14,15,18], 

Let  yh  denote  the  discrete  solution  space  for  the  velocity-pres- 
sure-level  set  triple  ju",  ph ,  (/)'’}  and  let  i^h  denote  the  discrete 
weighting  space  for  the  linear  momentum,  continuity  and  level 
set  equations  [w*,  qh,  t]h\.  The  RBVMS  formulation  of  the 
free-surface  equations  on  a  moving  domain  may  be  stated  as:  Find 
{u ",p*,  4 )*]  £  yh  such  that  V  (w*,  cj\  t/")  £  iVh\ 


'■'I 


(  <9u* 

W 


+  (u*  -  fi*)  ■  VX  -  f  )r/fi 


VAw"  :  er(u*,p*)c/fi  — 


w"  ■  c/h  1 


(r*)< 


X'Vv  ■  u'T/n 


tm  (  (u*  -  fi*)  ■  V,-w*  +  Va<? 


PT cVv  ■  w*rc(u",p*)dfi 

t 

TtfW*  ■  (] rM(u*,p" )  ■  V.vu *)c/fi 


■  Cm(u  h,ph)dQ. 


V,w" 


at  P 


:  (tmTm(u *,p"))  ®  (rMrM(u h,ph))dQ 


d£ 

W 


(u/!-fi *)■  V.vd»A  )c/£2 


(u*  -  fi*)  ■  V,-i/*  +  (u*  -  fi*)  ■  VX")  da  =  0 


(11) 


In  Eq.  (11)  all  the  integrals  are  taken  element-wise,  rM(u",  ph)  and 
rc{ u  ,  ph)  are  the  element-level  residuals  of  the  momentum  and 
continuity  equations  given  by  the  following: 


r*f(u  ,p)  =  p 


V.v-ff(u  ,p)  (12) 
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rc{»,p)  =  VA  ■  u 


(13) 

and  the  %'s  are  the  stabilization  parameters  discussed  in  the  sequel. 

In  the  presence  of  solid  walls,  (e.g.,  the  ship  hull),  it  is  benefi¬ 
cial  to  augment  the  coupled  weak  formulation  given  by  Eq.  (11) 
with  additional  terms  that  give  rise  to  weak  enforcement  of 
essential  boundary  conditions.  In  this  case,  the  left-hand  side  of 
Eq.  (1 1)  is  augmented  with  the  following  terms: 


satisfy  Eqs.  (20)-(21),  we  introduce  a  “pseudo-time”  variable  t  and 
integrate  the  following  semi-discrete  variational  equation  in 
pseudo-time:  Given  0;,(x,  0,  find  4>d  £  Vj ,  such  that  Vt/*  £  W* : 


n, 


^  +  a.vr^-sE(^)V/n 

r^a  ■  VAV  (jjf-  +  a  ■  V,0*  -  S.(0*)) 


r/£2  =  0 


(22) 


w''  ■  n(uh,pk)ndY 
(rg), 


where 


(2/tVXn  +  qh n)  ■  (u* 

(r«), 

w*-  ((u*-uA)  -n)(u* 

(I '8)7 

xBvth  ■  (u*  -  u*)t/r 

(rg), 

ug)dY 

Ug)dr 


(14) 


and 


a  =  S'£(0,')1 


(23) 


S£(0)=2//E(0)-l  (24) 


where  is  the  prescribed  fluid  velocity  on  (rg)(,  the  Dirichlet 
portion  of  the  fluid  domain  boundary,  and  (rg)j~  denotes  the 
“inflow”  part  of  (rg)f.  We  refer  the  readers  to  Refs.  [22,23,29]  for 
the  derivation  and  discussion  of  weak  boundary  conditions  for  the 
equations  of  fluid  mechanics.  %B  is  the  boundary  stabilization  or 
penalty  parameter. 

For  flows  at  very  high  Reynolds  number,  as  well  as  for  pure 
convection,  the  RBVMS  framework  alone  may  not  be  sufficiently 
robust.  In  this  case,  additional  discontinuity-capturing  terms  may 
be  added  to  the  left-hand  side  of  Eq.  (11) 


<  x  vy  *^ns  > 


VxrihKls\7xcl>hdn 


(15) 

(16) 


n, 


where  Kns  and  k/s  are  tensor-valued,  residual-based,  discontinuity¬ 
capturing  viscosities. 

The  density  and  viscosity  of  the  two-fluid  system  in  the  discrete 
setting  are  computed  as: 


^^(rt  +  P.fl-ftW*))  (17) 

p  =  pwH£(<t>h)  +  pa(l-H£(<t>h))  (18) 

where  H,.  is  the  regularized  Heaviside  function.  In  this  work,  we 
take: 

{0  if  0  <  — b 

^  (l +  y  +  ^sin(^))  if  |0|  <8  (19) 

1  if  0  >  £ 


are  an  effective  “convection”  velocity  and  a  regularized  sign  func¬ 
tion,  respectively. 

The  formulation  given  by  Eq.  (22)  is  the  SUPG  method  [30] 
applied  to  the  Eikonal  equation.  At  the  steady  state,  the  above 
problem  produces  a  new  level  set  field  0rf  with  the  signed  distance 
property  and  zero  level  set  coincident  with  that  of  0.  This  re¬ 
distancing  approach  was  first  proposed  in  Ref.  [4],  and  employed 
in  Refs.  [5,14,15]  for  finite  element  computations  of  free  surface 
phenomena. 

To  prevent  excessive  motions  of  the  air-water  interface  during 
the  re-distancing  procedure,  we  recommend  adding  a  penalty  term 
to  the  left-hand  side  of  Eq.  (22)  of  the  form: 


n, 


-  0 h)dSl 


(25) 


In  the  above  penalty  term,  the  constant  parameter  Xpen  is  multi¬ 
plied  by  //'(0j),  which  scales  as  l/h  in  the  interface  layer  render¬ 
ing  the  equations  dimensionally  consistent.  Furthermore,  //'(0j) 
is  identically  zero  away  from  the  interface;  hence  the  term  is  only 
active  where  it  is  necessary.  In  Sec.  5.1.1,  the  importance  of  this 
term  will  be  clearly  demonstrated. 

Both  convection  and  re-distancing  of  the  level  set  may  result  in 
the  loss  or  gain  of  the  total  fluid  mass.  The  amount  of  mass  deficit 
depends  on  many  factors.  One  is  more  likely  to  significantly  upset 
the  mass  balance  on  a  coarse  problem  mesh  than  on  a  fine  problem 
mesh.  Significant  mass  loss  or  gain  may  also  occur  when  the  dis¬ 
crete  equations  are  integrated  for  a  long  time  period.  In  this  case, 
seemingly  minor  mass  errors  for  a  given  time  step  may  compound 
into  a  large  mass  error  at  the  end  of  the  computation.  As  a  result, 
a  mass  correction  procedure  is  necessary.  We  begin  with  a  global 
mass  balance  law  for  a  moving  domain,  which  reads: 


where  e  ~  0(h)  defines  the  interface  width  between  the  air  and 
water  subdomains.  Note  that  £  — >  0  as  the  mesh  is  refined. 

Regularization  of  the  Heaviside  function  places  the  requirement 
on  the  level  set  function  07'  to  maintain  the  signed  distance  prop¬ 
erty.  This  means  that  the  value  of  <f>h  at  a  point  x  at  time  t  is  the 
perpendicular  distance  of  that  point  to  the  air-water  interface  Y™ . 
In  the  water  domain,  the  distance  takes  on  a  positive  value,  while 
in  the  air  domain  it  is  negative.  To  enforce  the  signed  distance 
property  of  the  level  set  function,  we  define  0^  such  that: 


d 

Jt 


n, 


pdCl  + 


u)  ■  nr/r  =  0 


(26) 


We  integrate  Eq.  (26)  in  time  over  the  time  step  ( tn ,  tn+i)  and  ap¬ 
proximate  the  term  corresponding  to  the  boundary  integral  using  a 
midpoint  rule  to  obtain: 


pn+1dQ 

&n+l 


p„dn 

Q„ 


Vjt02||  =  1  in  £2, 
05  =  0  on  C’ 


(20) 

(21) 


in  a  weak  sense.  Equation  (20)  is  the  Eikonal  partial  differential 
equation  subject  to  the  interior  constraint  given  by  Eq.  (21).  To 


+  Af„ 


Pn+l /2  (tt/j+1/2  —  fin+1/2)  ■  n„+i/2dY  —  0 

r„+ 1/2 


(27) 


which  is  a  time  discrete  version  of  the  mass  balance  equation.  To 
ensure  mass  balance  at  tn+i,  we  perturb  the  newly  re-distanced 
level  set  function  c pdh  by  a  global  constant  such  that  Eq.  (27)  is 
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satisfied  for  the  regularized  fluid  density  p  given  by  Eq.  (17).  This 
procedure  ensures  the  global  mass  balance  at  every  time  step.  To 
achieve  mass  balance  also  in  a  local  sense,  a  procedure  proposed 
in  Ref.  [14]  may  be  employed  instead. 

3.1.1  Computation  of  Mesh-Dependent  Parameters.  To 
define  the  method  parameters  that  depend  on  the  mesh  size,  we 
first  define  the  element  metric  tensor  as: 


used  to  regularize  the  fluid  density  and  viscosity  in  Eqs. 
(17) — ( 1 8).  The  structure  of  Eq.  (35)  gives  the  interface  half-width 
in  the  direction  of  Vx(j)h ,  which  is  orthogonal  to  the  interface,  as 
desired. 

Note  that  all  the  mesh-dependent  parameters  are  defined  in 
terms  of  the  element  metric  tensor  G,  which  automatically 
accounts  for  element  topology.  Alternative  definitions  of  stabiliza¬ 
tion  and  discontinuity-capturing  parameters  are  possible  and  may 
be  found  in  Refs.  [26,33,37-39], 


d <f  dj 

dx  dx 


(28) 


where  x  and  £  are  the  coordinates  of  the  physical  element  and  its 
parametric  counterpart,  respectively. 

The  RBVMS  formulation  parameters  zM,  zc  and  in  Eq.  (11) 
are  defined  as  (see  Ref.  [27] ) 


™  =  (^2  +  (u*  - fi/!)  ■  G(u"  -  + c>  (j)  G :  °r1/2  (29) 

TC  =  (ii-Gim)-1  (30) 

/  4  \  -1/2 

T^*(A?2  +  (U,!_fl^'G^_flAV  (31) 


3.2  Rigid  Body  Dynamics  Formulation.  We  begin  with  the 
definition  of  the  rigid  body  kinematics.  Let  Go  denote  the  refer¬ 
ence  configuration  of  the  rigid  object  and  let  G,  denote  its  current 
configuration.  Let  X  denote  the  reference  configuration  coordi¬ 
nates  of  the  rigid  body  and  let  x  denote  the  current  configuration 
coordinates  of  the  rigid  body.  The  rigid  body  kinematics  is 
composed  of  translation  and  rotational  motions  and  may  be 
summarized  as: 


x  —  Q(X  —  Xo)  +  Xo  +  do  (36) 

where  Q  is  a  rotation  matrix,  X0  is  a  center  of  mass,  and  do  is  the 
displacement  of  the  mass  center.  From  Eq.  (36)  we  obtain  the 
rigid  body  displacement  and  velocity  fields  as: 


where  C/  is  the  constant  emanating  from  the  element-wise  inverse 
estimate  (see  Ref.  [31]).  The  stabilization  parameter  z $  in 
Eq.  (22)  is  given  by: 


(4  VI/2 

T^=  (vAP  +  a'GaJ 


(32) 


and 


d  =  (Q  —  I)(X  —  X0)  +  d0  (37) 

d  =  Q(X  —  X0)  +  d0  (38) 


The  penalty  parameter  in  the  weak  Dirichlet  BC  formulation  zB  is 
given  by  (see  Ref.  [22]) 

zB  =  C^iVn  ■  Gn  (33) 


respectively. 

Let  x0  =  X0  +  do  denote  the  rigid  body  center  of  mass  in  the 
current  configuration.  Substituting  Eq.  (36)  in  Eq.  (38),  the  veloc¬ 
ity  field  may  be  expressed  in  terms  of  current  coordinates  as: 


where  C/  is  a  constant  emanating  from  the  boundary  inverse 
estimate.  An  alternative  definition  of  zB  that  is  based  on  the  wall 
function  ideas  may  be  found  in  Ref.  [23]. 

The  discontinuity  capturing  parameters  in  Eqs.  (15) — (16)  are 
given  by: 


t^ns  C  nx 


\\rM 


Kls  —  C  Is 


uref  VgTg 

&\ 

dt 


\M  +(uh-uh).\7xJ>h\ 


(34) 


V,r//  ■  GV,r// 


d  =  QQ  1  (x  -  x0)  +  d0  =  G(x  -  x0)  +  d0  (39) 


where 


G  =  QQ  1 


0  —  CO3  Ct>2 

co  3  0  — cui 

—CO  2  CO]  0 


(40) 


is  a  skew-symmetric  tensor  of  angular  velocities.  Note  that,  from 
Eq.  (40), 


Q  =  GQ 


(41) 


where  Cns  and  C/s  are  positive  constants  and  uref  is  a  reference  value 
of  flow  speed  (e.g.,  magnitude  of  inflow  velocity).  Note  that  the 
discontinuity-capturing  parameters  are  residual-based  and  isotropic. 
k,„  is  a  modified  version  of  the  YZ-/1  discontinuity  capturing  pro¬ 
posed  in  Refs.  [32,33]  corresponding  to  f  =  2,  which  is  the  least  in¬ 
trusive  definition.  The  definition  of  K/v  corresponds  to  the  CAU 
discontinuity  capturing  proposed  in  Ref.  [34],  and  may  also  be 
obtained  by  setting  /?  =  1  in  the  scalar  version  of  the  YZ-fi  disconti¬ 
nuity  capturing  method.  Alternative  definitions  of  discontinuity¬ 
capturing  terms,  which  are  not  residual  based,  can  be  found  in  Ref. 
[35]  in  the  context  of  incompressible  flows  and  in  Ref.  [36]  in  the 
context  of  incompressible  flows  with  thermal  coupling. 

Finally,  we  define  the  local  interface  half-width  e  as: 


which  is  a  key  relationship  that  we  will  employ  in  what  follows. 
Since  G  is  skew-symmetric,  it  has  the  axial  vector  co  given  by: 


V)\ 


0)2 


<03 


(42) 


and  the  rigid  body  velocity  field  given  by  Eq.  (39)  may  also  be 
written  as: 


d  =  co  x  (x  —  xq)  +  do 


(43) 


e  ma  "  "  (35) 

V  •  GV,4/’ 

where  a  is  an  integer  parameter  corresponding  to  the  number  of 
elements  used  across  the  air-water  interface.  This  definition  of  e  is 


In  3D,  the  translational  velocities,  do,  and  rotational  velocities,  co, 
are  the  six  degrees-of-freedom  that  define  the  kinematics  of  rigid 
body  motion. 

The  dynamics  of  the  rigid  body  is  governed  by  the  equations  of 
balance  of  global  linear  and  angular  momentum,  namely, 
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p&dSl  =  F 


d 

dt 


(x  —  x0)  x  pddQ  =  M 


(44) 

(45) 


where  F  and  M  are  the  global  force  and  moment  vectors  acting  on 
a  rigid  body  and  defined  by 


F  =  mg 


hr/r 


M  = 


(x  —  xq)  x  h  d r 


(46) 

(47) 


where  Ehm,  the  mesh  Young’s  modulus,  is  given  by: 

4  =  E"JJ  (55) 

Jx£  is  the  Jacobian  determinant  of  the  isoparametric  element  map¬ 
ping,  and  Em  and  vm  are  the  constant,  user-prescribed  nominal 
mesh  Young’s  modulus  and  Poisson  ratio.  We  typically  take 
Em  =  1  and  vm  =  0.3.  The  above  procedure  represents  the  so-called 
Jacobian  stiffening  method  (see  Ref.  [40]),  which  preserves  good 
mesh  quality  in  the  simulation.  In  particular,  the  small  elements, 
which  are  usually  placed  near  the  fluid-object  interface,  become 
“stiffer”  and  are  less  likely  to  deform  as  much  as  the  larger  ele¬ 
ments,  which  are  typically  placed  in  the  areas  where  the  solution 
is  not  expected  to  exhibit  complex  behavior. 


m  =  pdSl  is  the  mass  of  the  object,  g  is  the  gravity  vector,  and 
h  is  the  traction  vector  exerted  on  the  object  by  an  external  me¬ 
dium.  In  this  case,  h  is  the  fluid  traction  vector. 

Introducing  the  rigid  body  kinematics  into  the  linear  and  angu¬ 
lar  momentum  balance  equations,  we  obtain  the  following  system 
of  six  ordinary  differential  equations. 


^Mo)  =  F 

(48) 

^(J,m)  =  M 

(49) 

where  J,  is  the  inertia  tensor  in  the  current  configuration  given  by 

h  =  QJ0Qr  (50) 

with  J0  being  the  inertia  tensor  in  the  reference  configuration 
defined  as 


Jo  = 


P(X  -  Xo)  ■  (X 

Jn„ 

®  (X  -  X0)dSl 


x0)i  dn 


p(x  -  Xo) 

fio 


(51) 


Note  that  J0  is  a  function  of  the  reference  configuration  quantities 
only,  and  may  be  computed  once  and  for  all  for  a  given  rigid  object. 


4  Time  Integration  of  the  Free-Surface  Fluid-Object 
Interaction  Equations 

In  this  section  we  present  our  time  integration  algorithm  for  the 
coupled  free-surface  fluid-object  interaction  problem.  We  make 
use  of  the  generalized-a  time  integration  method  (see  Refs. 
[41,42])  for  the  free-surface  and  mesh  motion  equations.  The  rigid 
body  equations  are  integrated  with  a  midpoint  method.  The  latter 
choice  will  become  clear  in  the  sequel. 

The  generalize-a  method  applied  to  the  coupled  free-surface 
fluid-object  interaction  problem  may  be  recast  in  the  form  of  a 
three-stage  predictor-multicorector  algorithm  as  follows. 

Given  u,„  u„,  p„,  the  fluid  velocity,  its  time  derivative,  and  pressure 
nodal  degrees  of  freedom,  respectively,  <(),„  cfj,,,  the  level  set  and  its 
time  derivative  nodal  degrees  of  freedom,  respectively,  (do)„,  (do)„, 
Q,„  the  rigid  body  center  of  mass  displacement  and  velocity, 
angular  velocity,  and  the  rotation  matrix,  respectively,  and  d„,  u„,  it,, 
the  mesh  displacement,  velocity  and  acceleration  nodal  degrees  of 
freedom,  respectively,  at  time  level  /„,  find  the  conesponding  quanti¬ 
ties  at  time  level  4+t  by  executing  the  following  stages: 

Predictor  stage.  Initialize: 

t4+i,(o)  u  n 

y-  i . 

u«+l,(0)  —  — - — U« 


3.3  Motion  of  the  Fluid  Domain.  In  the  MITICT  framework, 
while  the  air-water  interface  is  captured  on  the  mesh  that  does  not 
conform  to  it,  the  fluid-rigid  object  interface  is  tracked  with  a  con¬ 
forming  mesh.  For  this,  the  computed  rigid  object  displacement 
and  velocity  (see  Eqs.  (37)  and  (38))  are  imposed  as  boundary  con¬ 
ditions  for  the  motion  of  the  fluid  domain  mesh.  The  fluid  mesh  dis¬ 
placement  in  the  interior  of  the  fluid  domain  is  found  using  the 
equations  of  linear  elastostatics:  Find  d;'(0  such  that  Vw* 


*wVvj(d \t)  -  d*(/)) 


d  Q 


Vi  •  w' 


'khVi  ■  (d\t)  -  d*(f))dQ  =  0  (52) 


-  h 

In  Eq.  (52)  d  (f)  and  £2,  are  the  mesh  displacement  and  mesh  con¬ 
figuration  at  t  <  t,  both  considered  known.  In  practice,  we  take  Si; 
to  be  the  configuration  at  the  previous  time  step.  This  renders  the 
mesh  motion  problem  linear  in  the  unknown  displacement  vari¬ 
able  d/7(0-  Also  in  Eq.  (52)  //'  and  kh  are  the  mesh  Lame  parame¬ 
ters  chosen  to  be: 


p>  = 


2(1 


vmEi 


(1  +  vm)(l  -  2vm) 


(53) 

(54) 


Pn+1,(0)  P  n 
<t>n+l,(0)  =  4*n 

•k+lJO)  =  “  4*n 

U„+1,(0)  =  fin 

-  y-  1  - 

u«+l,(0)  —  — - — U« 

At2 

d„+i,(o)  =  d„  4-  A?u„  +  ~2~  ((1  —  2/l)u„  +  2/Ju„+ij(o)) 

The  subscript  0  on  the  left-hand  side  quantities  is  the  iteration 
index,  which  is  set  to  0  at  the  predictor  stage. 

Multicorrector  stage.  Repeat  the  following  iterations  for 
/=  1,2 ...lmax,  where  I  is  the  iteration  index  and  lmax  is  the  maximum 
allowable  number  of  nonlinear  iterations  set  for  this  time  step: 

1.  (a)  Evaluate  iterates  at  the  intermediate  time  levels  as: 

t4+a/,(/)  tln  +  rJ-j  (utt  i  |  ('/  i  j  U„) 

■  ’J-f.U)  —  On  +  rJ-,n  (U/j  |  | .(/  |  j  U // ) 

<t>«+ar,(/)  =  *}>/!  +  “/■(<I>«+1,(/-1)  “  <t>/7) 

(56) 

<J)n+a„,(/)  =  <t)„  +  «m(<IV|-l,(Z-l)  “  (t>„) 

d/7+tz/.,(/)  d^  — F  y-f  ( d/j  i  i  y  i  ^  d/;) 
d/7+ct/.,(/)  d/j  — F  y.m  (d„+i,(/_i)  -  d„) 
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(. b )  Use  the  intermediate  solutions  to  assemble  R'"0"',  Rcon,  and 
Rk,  the  discrete  residuals  of  the  momentum,  continuity  and 
level-set  equations,  and  the  corresponding  matrices  of  the  linear 
equation  system, 


dRmom  .  .  dR"wm 

- Au  -i - 

dUn+l 


dRcon 

9u„+i 

dRls 

■n+ 1 


dir 


Au 


Au 


aR  mom 

Ap  H - t - A(j>  = 

dp„+1  d§n\\ 


R 


dRcon  ,  dRcon 

- Ap  ■+- - 

dp„+i 

dRh 

W„ 


d< t>„ 

A<i> = -r: 


A<j>  =  -r; 


(0 


(0 


Solve  this  linear  system  using  a  preconditioned  GMRES  algorithm 
(see  Ref.  [43])  to  a  specified  tolerance.  We  note  that  the  left-hand 
side  matrix  couples  the  fluid  and  level  set  degrees  of  freedom. 
This  strong  coupling  between  the  fluid  and  level  set  equations  allows 
us  to  use  larger  time  step  sizes  without  the  risk  of  encountering  an 
instability  associated  with  staggered  approaches  in  which  the  fluid 
solve  is  followed  by  a  level  set  solve. 

(c)  Update  the  fluid  and  level  set  solution  as: 

>in+i,(/)  =  u„+i,{;)  +  Au 

u„+i,(/)  =  +  yA?Au 

P„+i,(t)  =  P«+t,(i)  +  AP  (57) 

<k+t,(i)  =  <k+i,(/)  +  Acj> 

4,n+l,(/)  =  4*n+l,(/)  +  yAfA«j> 


(. b )  Use  the  half-step  fluid  solution  to  assemble  the  F"+1/2  and 
M'!+1/2,  force  and  moment  load  on  the  rigid  object,  respectively, 
using  Eqs.  (46)  and  (47).  Using  these  loads,  we  solve  the  follow¬ 
ing  discrete  rigid  body  equations  for  (d0), ,+],(/)>  (do)„+1  m, 
•U/i+t,©!  and  Qn+i.qy 


(do) 


n+l,(()  (UOln 

At 


=  F"+5 


(do)«+i,(/)  -  do)„  1  /  ■  ■ 

•  2  y(d°)«+t,(/)  +  (do)„ 


At 


Jfl+l,(/)m«+l,(f)  J /i<a«  _  jyjn-fi 


At 


'  “  -  4  (^n+l.W  +  (Q»+1  ,(/)  +  Qh) 


(59) 

(60) 
(61) 
(62) 


where  Eq.  (50)  is  used  to  evaluate  the  inertia  tensor  as: 

=  Q„+i,(,)JoQl+1,w  (63) 


Note  that  Eqs.  (61) — (63)  are  nonlinear  and  coupled,  requiring  iter¬ 
ation.  We  solve  Eq.  (61)  first,  where  we  lag  the  inertia  tensor.  We 
solve  Eq.  (62)  next  and  use  the  new  rotation  matrix  to  update  the 
inertia  tensor  in  Eq.  (63).  This  procedure  gives  convergence  in  a 
few  iterations.  The  readers  should  keep  in  mind  that  these  iteration 
are  very  inexpensive. 

3.  (a)  Given  the  displacement  of  the  rigid  object,  use  it  to  set 
the  boundary  conditions  for  the  fluid  mesh  motion  problem  as: 


2  (a)  Evaluate  the  updated  fluid  solution  at  the  half-step  as: 


d, ,+!,(/)  -  (Q„+i,(/)  -  l)  (x  -  X0)  +  (do)„+ii(,)  (64) 


u«+i  (/)  =  j  (Un  +  u«+i.w) 

“„+!(/)  =  j  (“«  + 

<hi+l(()  =  2  (‘t'n  +  <t,n+l,(/)) 

‘i’n+l.to  =  2  (<j>»  +  <Ki+l,(o) 

d„+i  (/)  =  j  (d«  +  d„+i,(/)) 

u„+i,(z)  =  2  (u«  +  u„+i,(/)) 


(h)  Evaluate  the  intermediate  time  level  mesh  displacement  as: 

d„  i  y.f.(j)  d„  T  y~f(dn  \  \ij  i  j  d„)  (65) 

(c)  Use  the  intermediate  mesh  displacement  solution  to  assem- 
(58)  Rmesh,  the  discrete  residual  of  the  mesh  motion  problem, 

and  the  corresponding  left-hand  side  matrix  of  the  linear  equa¬ 
tion  system: 


rrnsh 

<9u„+i 


Au  = 


o  mesh 
“K(0 


(66) 


Fig.  2  Dam  break  with  obstacle.  Problem  setup. 
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Fig.  3  Dam  break  with  obstacle.  Snapshots  of  the  water  subdomain  colored  by 
the  fluid  speed  at  f=2.0s.  Top:  solution  without  penalty.  Bottom:  solution  with 
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Apen  1  ■ 


Solve  this  linear  system  using  a  preconditioned  Conjugate 
Gradient  algorithm  to  a  specified  tolerance. 

(d)  Update  the  mesh  motion  solution  as: 

Un+t,(t)  =  tin+t,(/)  +  Au 

u«+i,(0  =  uB+t,(t)  +  yArAu  (67) 

=  d„+i,(;)  +  /?At2Ah 

Level  Set  Correction  Stage.  Re-distance  the  level  set  field  accord¬ 
ing  to  Eq.  (22)  and  correct  for  mass  deficit  according  to  Eq.  (27). 

This  completes  the  time  step,  at  which  point  the  time  step  coun¬ 
ter  is  incremented  and  we  go  back  to  the  predictor  stage. 

Remark.  Our  rigid  body  time  integration  algorithm  amounts 
to  applying  the  midpoint  rule  to  Eqs.  (48)  and  (49)  as  well  as 
Eq.  (41).  That  is,  the  rotation  matrix  Q  is  not  explicitly  computed 
from  the  rotation  angles,  but  is  carried  as  a  separate  problem 
unknown.  The  choice  of  the  midpoint  time  integration  for  Eq.  (41) 
results  in  the  following  remarkable  property:  if  Q^Q„ 
=  Q„Ql  =  I  then  Q,'  , Q,, ,  i  =  Q„ .  |Q,',  [  =  I.  That  is,  once  ini¬ 


tialized  as  a  proper  rotation  matrix,  Q  remains  a  proper  rotation 
matrix  during  the  entire  simulation.  This  result  is  due  to  Ref.  [44]. 

5  Numerical  Results 

In  this  section,  we  present  two  numerical  examples.  The  first 
one  is  the  well-known  MARIN  dam  break  problem  [14—17], 
which  is  a  free-surface  flow  without  fluid-object  interaction.  The 
problem  is  solved  on  a  fixed  mesh  and  is  used  to  illustrate  the 
effect  of  the  penalty  term  in  the  re-distancing  step  given  by 
Eq.  (25).  This  problem  is  also  used  to  compare  a  fully-coupled 
versus  staggered  approach  to  level  set  convection.  The  second 
example  makes  use  of  a  DTMB  5415  Navy  combatant  at  lab  scale 
in  head  waves  of  large  amplitude. 

5.1  Dam  Break  With  Obstacle.  The  problem  consists  of  a 
1.22m  x  1  m  x  0.55  m  column  of  water,  initially  at  rest,  that 
collapses  under  the  action  of  gravity  and  impacts  a  stationary 
object.  The  computational  domain  is  a  rectangular  box  with 
dimensions  3.22m  xlmxlm.  The  object  has  dimensions 


Fig.  4  Dam  break  with  obstacle.  Time  series  of  the  pressure  at  four  locations  on 
the  obstacle  (see  Fig.  2).  Very  close  correlation  with  experimental  data  is  obtained. 
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0.161  m  x  0.161  m  x  0.403  m  and  is  placed  at  the  back  end  of  the 
tank.  The  problem  setup  is  shown  in  Fig.  2.  Experiments  for  this 
test  case  were  performed  at  the  Maritime  Research  Institute 
Netherlands  (MARIN),  and  the  data  is  often  used  to  validate  free- 
surface  software  for  marine  engineering  applications. 

In  Refs.  [14,15]  we  computed  this  test  case  with  linear  hexahe- 
dral  and  quadratic  NURBS  and  reported  very  good  correlation 
between  experimental  and  computational  results.  In  this  paper,  we 
employ  linear  tetrahedral  elements  and  use  the  MARIN  example 
to  assess  the  influence  of  the  penalty  parameter  in  the  re¬ 
distancing  equations,  and  illustrate  the  increased  robustness  of  the 
direct  coupling  between  the  fluid  flow  and  level  set  equations. 

5.1.1  Influence  of  Penalty  in  the  Re-Distancing.  Simulations 
with  two  different  penalty  parameters,  7pE„  =  0,  ipE„  =  1,  are  per¬ 
formed  on  a  mesh  consisting  of  344,401  tetrahedral  elements  and 
60,797  nodes.  The  simulation  was  run  until  T  =6.0.t  with  a  time 
step  of  At  =  0.01.  Snapshots  of  the  solution  for  both  simulations 
are  shown  in  Fig.  3.  At  f  =  2.0.t,  the  water  has  already  hit  the  ob¬ 
stacle  and  rebounded  off  the  back  wall  of  the  tank  to  form  a  break¬ 
ing  wave.  The  case  with  penalty  shows  a  crisp  resolution  of  the 
air-water  interface  and  is  correctly  predicting  the  formation  of  a 
breaking  wave.  The  no-penalty  case,  however,  produced  an  over¬ 
diffuse  interface  and;  as  a  result,  it  is  unable  to  predict  wave 
breaking.  This  example  clearly  illustrates  the  fact  that  the  regular¬ 
ized  sign  function  given  by  Eq.  (24)  is  not  a  sufficient  mechanism 
by  itself  to  keep  the  interface  intact  during  re-distancing. 

Figure  4  shows  the  pressure  history  on  four  different  locations 
on  the  obstacle.  The  simulated  results  correlate  quite  well  with  the 
experiments  for  such  a  coarse  discretization. 


5.1.2  Fully-Coupled  Versus  Staggered  Fluid-Level  Set 
Simulations.  Here  we  solve  a  MARIN  dam  break  problem  using 
the  solution  strategy  presented  in  Sec.  4  and  compare  it  with  the 
solution  strategy  we  presented  in  Refs.  [14,15].  The  latter  decou¬ 
ples  the  Navier-Stokes  and  level  set  convection  equations  and  per¬ 
forms  one  convection  step,  which  includes  re-distancing  and  mass 
correction,  at  the  end  of  the  time  step.  This  methodology  is  com¬ 
monly  used  in  free  surface  computations  reported  in  the  literature. 
For  this  test,  we  reduce  the  mesh  size  to  103  251  tetrahedral  ele¬ 
ments  and  19  071  nodes.  We  increase  the  time  step  size  to 
At  =0.025  and  compute  the  problem  until  T=  12.5  5.  By  this 
time,  the  fluid  solution  should  be  well  on  its  way  to  a  steady  result. 
Figure  5  shows  the  fluid  solution  at  t=  12.5.?,  the  final  time.  The 
fully  coupled  strategy  gives  a  physical,  near  steady-state  result, 
while  the  staggered  case  predicts  large-magnitude  sloshing,  which 
is  unphysical.  This  is  due  to  the  large  time  step  size  that  triggers 
an  instability  in  the  simulations.  This  example  illustrates  the  dan¬ 
ger  of  using  staggered  methods  with  time  step  sizes  that  exceed 


Fig.  5  Dam  break  with  obstacle.  Snapshots  of  the  water  sub- 
domain  colored  by  the  fluid  speed  at  f=  12.5s.  Top:  fully 
coupled  simulation.  Bottom:  staggered  simulation.  For  the  cho¬ 
sen  time  step  size,  the  coupled  simulation  produces  a  physical, 
near  steady-state  result,  while  the  staggered  approach  gives 
unphysical,  large-magnitude  sloshing. 


the  stability  limits  of  the  formulation  and  lead  to  seemingly  con¬ 
vergent,  yet  completely  unphysical  results. 


5.2  Ship  in  Head  Sea.  Here  we  simulate  the  DTMB  5415 
Navy  combatant  at  lab  scale.  This  ship  has  been  investigated  by 
other  researchers,  both  experimentally  and  computationally  (see 
Refs.  [45-47]).  The  length  of  the  ship  hull  is  5.12  m.  The  ship 
mass,  center  of  gravity  and  inertia  tensor  are  computed  by  mesh¬ 
ing  the  ship  interior  and  performing  a  direct  computation.  The 
total  ship  volume  is  1.366m3.  The  ship  mass  is  equal  to  532.3  kg. 
It  is  obtained  by  multiplying  the  volume  of  the  ship  below  the 
water  line  by  the  constant  water  density.  The  center  of  gravity  and 
the  inertia  tensor  are  computed  assuming  the  ship’s  effective  den¬ 
sity  (i.e.,  the  ship  mass  divided  by  its  total  volume),  which  results 
in 
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b)  t=9.50s 


Fig.  6  DTMB  5415  in  head  sea.  Snapshots  of  the  ship  negotiat¬ 
ing  high-amplitude  waves.  The  water  surface  is  colored  by  the 
fluid  speed. 
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Fig.  7  DTMB  5415  in  head  sea.  Time  history  of  ship  motion. 


(a)  Force  load  on  ship  hull 


(b)  Moment  load  on  ship  hull 

Fig.  8  DTMB  5415  in  head  sea.  Time  history  of  forces  and  moments  acting  on  the  ship. 
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We  compute  the  ship  in  head  waves,  meaning  the  waves  that 
travel  in  the  direction  opposite  to  that  of  the  ship.  We  assume  that 
the  ship  speed  is  Uin  =  1.873  mis,  which  gives  the  Froude  number 
of  0.25  based  on  the  ship  length. 

We  make  use  of  the  linear  Airy  waves  [48]  to  prescribe  inlet 
boundary  conditions.  The  Airy  waves  may  be  derived  using  poten¬ 
tial  theory,  and  are  specified  as  follows:  Given,  the  wave  ampli¬ 
tude,  wave  length  and  water  depth,  Aw  =  0.2m,  L„.  =  5.72  m  and 
/;  =  3.49  m,  respectively,  we  compute  k  =  2n/Lw,  the  angular 
wavenumber,  m  =  gktanh(kh) :  the  wave  phase  speed,  and 
Av  =  the  velocity  amplitude.  With  these  definitions,  the 

Airy  waves  are  given  by: 


u  =  Av  cosh(fe)  cos(fcv  —  mi)  +  Uj„ 

(70) 

v  =  0 

(71) 

w  =  Av  sinh(Lz)  sin(fct  —  cot ) 

(72) 

< j>  =  Am  cos(fct  —  mt)  +  h  —  z 

(73) 

where  u  =  (u,v,w)  is  the  fluid  velocity  vector  and  the  air-water 
interface  in  the  hydrostatic  configuration  is  assumed  to  be  located 
at  z  =  0. 

In  the  continuum  setting,  the  fluid  traction  vector  h  is  defined 
as: 

h  =  pa  —  2/rVvu  ■  n  (74) 

In  the  discrete  setting,  in  the  presence  of  weakly-enforced  bound¬ 
ary  conditions,  basic  conservation  arguments  (see  Ref.[49])  lead 
to  the  following  definition  of  the  traction  vector 

h  =  ph n  —  2/jV'u*  ■  n  +  zB( nh  —  d)  (75) 

We  use  this  definition  of  h  for  the  computation  of  global  force  and 
moment  vectors  (see  Eqs.  (46)  and  (47))  acting  on  the  ship  hull. 

The  simulation  was  performed  on  a  mesh  consisting  of 
6  285  445  linear  tetrahedral  elements  and  1  059  174  nodes.  The 
simulation  took  5000  time  steps  at  a  fixed  time  step  size  of 
Af=  0.0025  s.  The  ship  was  allowed  to  move  vertically,  to  pitch 
and  to  roll,  while  the  rest  of  the  rigid  body  DOFs  were 
constrained. 

Figure  6  shows  the  snapshots  of  the  ship  negotiating  high- 
amplitude  waves.  The  bottom  part  of  Fig.  6  shows  the  ship  par¬ 
tially  submerged  in  water,  which  is  a  result  of  the  oncoming  wave 
hitting  the  bow  of  the  ship.  In  this  case,  near  the  bow,  the  free  sur¬ 
face  experiences  topological  changes,  which  necessitates  the  use 
of  an  interface-capturing  method  for  this  class  of  problems. 

The  time  history  of  the  ship  motion  is  given  in  Fig.  7.  Although 
the  prescribed  Airy  waves  only  have  a  single  frequency,  multiple 
frequencies  are  present  in  the  ship’s  response.  Note  that  the  ship 
develops  a  low-amplitude,  chaotic  rolling  motion. 

Figure  8(a)  shows  the  time  history  of  the  thrust  force  necessary 
to  maintain  the  ship  moving  forward  at  constant  speed.  The  time 
history  of  the  forces  and  moments  in  the  unconstrained  directions 
are  shown  in  Figs.  8 (h)-(d). 

6  Conclusions 

A  free  surface-rigid-body  interaction  modeling  approach  is  pre¬ 
sented  using  the  MITICT  approach,  in  which  the  air-water  inter¬ 
face  is  handled  by  means  of  an  interface-capturing  level  set 
method,  and  the  fluid-rigid  body  interface  is  handled  with  an 
interface-tracking  ALE  method.  The  RBVMS  formulation  is  used 
for  the  fluid  mechanics  and  level  set  equations.  Weak  enforcement 
of  essential  boundary  conditions  is  employed  on  all  no-slip 
surfaces  for  better  approximation  of  thin  fluid  boundary  layers. 
The  rigid-body  time  integration  algorithm  is  proposed  where  a 
separate  time  evolution  equation  is  solved  for  the  rotation  matrix, 
which  becomes  an  additional  problem  unknown.  The  use  of  the 


midpoint  time  integration  algorithm  for  the  rotation  matrix  results 
in  the  exact  discrete  preservation  of  its  orthonormal  property,  a 
result  that  is  due  to  Ref.  [44].  In  the  proposed  methodology,  the 
strong  coupling  between  the  fluid  and  level  set  equations  at  the 
Newton  iteration  level  allows  us  to  robustly  march  in  time  by 
using  larger  time-steps  than  in  the  more  usual  staggered  (or  split 
operator)  methods.  We  also  illustrate  the  importance  of  the  pen¬ 
alty  parameter  in  the  level  set  re-distancing  equations  for  the  pres¬ 
ervation  of  the  sharpness  of  the  air-water  interface. 

In  future  work,  we  plan  to  enhance  the  structural  modeling  in 
our  framework  by  also  considering  elastic  bodies.  We  likewise 
plan  to  simulate  ships  at  larger  spatial  scales,  and  look  at  other 
applications,  such  as  modeling  of  offshore  wind  turbines. 


Acknowledgment 

This  research  was  supported  in  part  by  an  appointment  to  the 
Postgraduate  Research  Participation  Program  at  the  U.S.  Army 
Engineering  Research  and  Development  Center,  Coastal  and 
Hydraulics  Laboratory  (ERDC-CHL)  administered  by  the  Oak 
Ridge  Institute  for  Science  and  Education  through  an  interagency 
agreement  between  the  U.S.  Department  of  Energy  and  ERDC. 
The  corresponding  author  was  supported  through  ARO  Award 
W911NF-11-1-0083.  Permission  was  granted  by  the  Chief  of 
Engineers  to  publish  this  information. 

References 

[1]  Takizawa,  K.,  Tanizawa,  K.,  Yabe,  T.,  and  Tezduyar,  T.,  2007,  “Ship  Hydrody¬ 
namics  Computations  With  the  CIP  Method  Based  on  Adaptive  Soroban 
Grids,”  Int.  J.  Numer.  Methods  Fluids,  54,  pp.  101 1-1019. 

[2]  Tezduyar,  T.,  Aliabadi,  S.,  and  Behr,  M.,  1998,  “Enhanced-Discretization 
Interface-Capturing  Technique  (EDICT)  for  Computation  of  Unsteady  Flows 
With  Interfaces,”  Comput.  Methods  Appl.  Mech.  Eng.,  155,  pp.  235-248. 

[3]  Sethian,  J.,  1999,  Level  Set  Methods  and  Fast  Marching  Methods,  Cambridge 
University  Press,  London. 

[4J  Sussman,  M.,  Smereka,  P.,  and  Osher,  S.,  1994,  “A  Level  Set  Approach  for 
Computing  Solutions  to  Incompressible  Two-Phase  Flows,”  J.  Comput.  Phys., 
114,  pp.  146-159. 

[5]  Nagrath,  S.,  Jansen,  K.,  and  Lahey,  R.,  2005,  “Computation  of  Incompressible 
Bubble  Dynamics  With  a  Stabilized  Finite  Element  Level  Set  Method,”  Com¬ 
put.  Methods  Appl.  Mech.  Eng.,  194,  pp.  4565^-587. 

[6]  Olsson,  E.,  and  Kreiss,  G.,  2005,  “A  Conservative  Level  Set  Method  for  Two 
Phase  Flow,”  J.  Comput.  Phys.,  210,  pp.  225-246. 

[7]  Tezduyar,  T.,  2001, “Finite  Element  Methods  for  Flow  Problems  With  Moving 
Boundaries  and  Interfaces,”  Arch.  Comput.  Methods  Eng.,  8,  pp.  83-130. 

[8]  Akin,  J.,  Tezduyar,  T.,  and  Ungor,  M.,  2007,  “Computation  of  Flow  Problems 
With  the  Mixed  Interface-Tracking/Interface-Capturing  Technique  (MITICT),” 
Comput.  Fluids,  36,  pp.  2—11. 

[9]  Tezduyar,  T.,  Behr,  M.,  and  Liou,  J.,  1992,  “A  New  Strategy  for  Finite  Element 
Computations  Involving  Moving  Boundaries  and  Interfaces — The  Deforming- 
Spatial-Domain/Space-Time  Procedure:  I.  The  Concept  and  the  Preliminary 
Numerical  Tests,”  Comput.  Methods  Appl.  Mech.  Eng.,  94(3),  pp.  339-351. 

[10]  Tezduyar,  T.,  Behr,  M.,  Mittal,  S.,  and  Liou,  J.,  1992,  “A  New  Strategy  for  Fi¬ 
nite  Element  Computations  Involving  Moving  Boundaries  and  Interfaces — The 
Deforming-Spatial-Domain/Space-Time  Procedure:  II.  Computation  of  Free- 
Surface  Flows,  Two-Liquid  Flows,  and  Flows  With  Drifting  Cylinders,”  Com¬ 
put.  Methods  Appl.  Mech.  Eng.,  94(3),  pp.  353-371. 

[11]  Cruchaga,  M.,  Celentano,  D.,  and  Tezduyar,  T.,  2007,  “A  Numerical  Model 
Based  on  the  Mixed  Interface-Tracking/Interface-Capturing  Technique 
(MITICT)  for  Flows  With  Fluid-Solid  and  Fluid-Fluid  Interfaces,”  Int.  J. 
Numer.  Methods  Fluids,  54,  pp.  1021-1031. 

[12]  Cruchaga,  M.,  Celentano,  D.,  and  Tezduyar,  T.,  2001,  “A  Moving  Lagrangian 
Interface  Technique  for  Flow  Computations  Over  Fixed  Meshes,”  Comput. 
Methods  Appl.  Mech.  Eng.,  191,  pp.  525-543. 

[13]  Hughes,  T.  J.  R.,  Liu,  W.  K.,  and  Zimmermann,  T.  K.,  1981,  “Lagrangian- 
Eulerian  Finite  Element  Formulation  for  Incompressible  Viscous  Flows,”  Com¬ 
put.  Methods  Appl.  Mech.  Eng.,  29,  pp.  329-349. 

[14]  Kees,  C.,  Akkerman,  I.,  Farthing,  M.,  and  Bazilevs,  Y.,  2011,  “A  Conservative 
Level  Set  Method  Suitable  for  Variable-Order  Approximations  and  Unstruc¬ 
tured  Meshes,”  J.  Comput.  Phys.,  230,  pp.  4536-4558. 

[15]  Akkerman,  I.,  Bazilevs,  Y.,  Kees,  C.,  and  Farthing,  M.,  2011,  “Isogeometric 
Analysis  of  Free-Surface  Flow,”  J.  Comput.  Phys.,  230,  pp.  4137-^4152. 

[16]  Kleefsman,  K.,  Fekken,  G.,  Veldman,  A.,  Iwanowski,  B.,  and  Buchner,  B., 
2005,  “A  Volume-of-Fluid  Based  Simulation  Method  for  Wave  Impact  Prob¬ 
lems,”  J.  Comput.  Phys.,  206,  pp.  363-393. 

[17]  Elias,  R.,  and  Coutinho,  A.,  2007,  “Stabilized  Edge-Based  Finite  Element  Sim¬ 
ulation  of  Free-Surface  Flows,”  Int.  J.  Numer.  Methods  Fluids,  54,  pp. 
965-993. 


010905-10  /  Vol.  79,  JANUARY  2012 


Transactions  of  the  ASME 


[18]  Lins,  E.  F.,  Elias,  R.  N.,  Rochinha,  F.  A.,  and  Coutinho,  A.  L.  G.  A.,  2010, 
“Residual-Based  Variational  Multiscale  Simulation  of  Free  Surface  Flows,” 
Comput.  Mech.,  46,  pp.  545-557. 

[19]  Bazilevs,  Y.,  Calo,  V.,  Cottrel,  J.,  Hughes,  T.  J.  R.,  Reali,  A.,  and  Scovazzi,  G., 
2007,  “Variational  Multiscale  Residual-Based  Turbulence  Modeling  for  Large 
Eddy  Simulation  of  Incompressible  Flows,”  Comput.  Methods  Appl.  Mech. 
Eng.,  197,  pp.  173-201. 

[20]  Bazilevs,  Y.,  Calo,  V.,  Zhang,  Y.,  and  Hughes,  T.  J.  R.,  2006,  “Isogeometric 
Fluid-Structure  Interaction  Analysis  With  Applications  to  Arterial  Blood 
Flow,”  Comput.  Mech.,  38,  pp.  310-322. 

[21]  Akkerman,  I.,  Bazilevs,  Y.,  Calo,  V.,  Hughes,  T.  J.  R.,  and  Hulshoff,  S.,  2008, 
“The  Role  of  Continuity  in  Residual-Based  Variational  Multiscale  Modeling  of 
Turbulence,”  Comput.  Mech.,  41,  pp.  371-378. 

[22]  Bazilevs,  Y.,  Michler,  C.,  Calo,  V.,  and  Hughes,  T.  J.  R.,  2007,  “Weak  Dirichlet 
Boundary  Conditions  for  Wall-Bounded  Turbulent  Flows,”  Comput.  Methods 
Appl.  Mech.  Eng.,  196,  pp.  4853—4862. 

[23]  Bazilevs,  Y.,  Michler,  C.,  Calo,  V.,  and  Hughes,  T.  J.  R.,  2010,  “Isogeometric 
Variational  Multiscale  Modeling  of  Wall-Bounded  Turbulent  Flows  With 
Weakly-Enforced  Boundary  Conditions  on  Unstretched  Meshes,”  Comput. 
Methods  Appl.  Mech.  Eng.,  199(13-16),  pp.  780-790. 

[24]  Bazilevs,  Y.,  Gohean,  J.,  Hughes,  T.  J.  R.,  Moser,  R.,  and  Zhang,  Y.,  2009, 
“Patient-Specific  Isogeometric  Fluid-Structure  Interaction  Analysis  of  Thoracic 
Aortic  Blood  Flow  Due  to  Implantation  of  the  Jarvik  2000  Left  Ventricular 
Assist  Device,”  Comput.  Methods  Appl.  Mech.  Eng.,  198,  pp.  3534-3550. 

[25]  Bazilevs,  Y.,  Calo,  V.,  Hughes,  T.  J.  R.,  and  Zhang,  Y.,  2008,  “Isogeometric 
Fluid-Structure  Interaction:  Theory,  Algorithms,  and  Computations,”  Comput. 
Mech.,  43,  pp.  3—37. 

[26]  Hsu,  M.,  Bazilevs,  Y.,  Calo,  V.,  Tezduyar,  T.,  and  Hughes,  T.  J.  R.,  2010, 
“Improving  Stability  of  Multiscale  Formulations  of  Fluid  Flow  at  Small  Time 
Steps,”  Comput.  Methods  Appl.  Mech.  Eng.,  199,  pp.  828-840. 

[27]  Bazilevs,  Y.,  Hsu,  M.-C.,  Akkerman,  I.,  Wright,  S.,  Takizawa,  K.,  Henicke,  B., 
Spielman,  T.,  and  Tezduyar,  T.,  2011,  “3D  Simulation  of  Wind  Turbine  Rotors 
at  Full  Scale.  Part  I:  Geometry  Modeling  and  Aerodynamics,”  Int.  J.  Numer. 
Methods  Fluids,  65,  pp.  207-235. 

[28]  Bazilevs,  Y.,  Hsu,  M.-C.,  Kiendl,  J.,  Wuechner,  R.,  and  Bletzinger,  K.-U., 
2011,  “3D  Simulation  of  Wind  Turbine  Rotors  at  Full  Scale.  Part  II:  Fluid- 
Structure  Interaction,”  Int.  J.  Numer.  Methods  Fluids,  65,  pp.  236—253. 

[29]  Bazilevs,  Y.,  and  Hughes,  T.  J.  R.,  2007,  “Weak  Imposition  of  Dirichlet  Bound¬ 
ary  Conditions  in  Fluid  Mechanics,”  Comput.  Fluids,  36,  pp.  12-26. 

[30]  Brooks,  A.,  and  Hughes,  T.  J.  R.,  1982,  “Streamline  Upwind/Petrov-Galerkin 
Formulations  for  Convection  Dominated  Flows  With  Particular  Emphasis  on 
the  Incompressible  Navier-Stokes  Equations,”  Comput.  Methods  Appl.  Mech. 
Eng.,  32,  pp.  199-259. 

[31]  Harari,  I.,  and  Hughes,  T.  J.  R.,  1992,  “What  are  C  and  h?:  Inequalities  for  the 
Analysis  and  Design  of  Finite  Element  Methods,”  Comput.  Methods  Appl. 
Mech.  Eng.,  97,  pp.  157-192. 

[32]  Tezduyar,  T.,  2004,  “Finite  Element  Methods  for  Fluid  Dynamics  With  Moving 
Boundaries  and  Interfaces,”  Encyclopedia  of  Computational  Mechanics,  E. 
Stein,  R.  D.  Borst,  and  T.  J.  R.  Hughes,  eds.,  Vol.  3:  Fluids,  John  Wiley  & 
Sons,  New  York. 


[33]  Tezduyar,  T.,  and  Senga,  M.,  2006,  “Stabilization  and  Shock-Capturing  Param¬ 
eters  in  SUPG  Formulation  of  Compressible  Flows,”  Comput.  Methods  Appl. 
Mech.  Eng.,  195,  pp.  1621—1632. 

[34]  Galeao,  A.,  and  do  Carmo,  E.  D.,  1988,  “A  Consistent  Approximate  Upwind 
Petrov-Galerkin  Method  for  Convection-Dominated  Problems,”  Comput. 
Methods  Appl.  Mech.  Eng.,  68(1),  pp.  83-95. 

[35]  Rispoli,  F.,  Corsini,  A.,  and  Tezduyar,  T.,  2007,  “Finite  Element  Computation 
of  Turbulent  Flows  With  the  Discontinuity-Capturing  Directional  Dissipation 
(DCDD),”  Comput.  Fluids,  36,  pp.  121-126. 

[36]  Tezduyar,  T.,  Ramakrishnan,  S.,  and  Sathe,  S.,  2008,  “Stabilized  Formulations 
for  Incompressible  Flows  With  Thermal  Coupling,”  Int.  J.  Numer.  Methods 
Fluids,  57,  pp.  1189-1209. 

[37]  Tezduyar,  T.,  and  Osawa,  Y.,  2000,  “Finite  Element  Stabilization  Parameters 
Computed  From  Element  Matrices  and  Vectors,”  Comput.  Methods  Appl. 
Mech.  Eng.,  190,  pp.  41 1^130. 

[38]  Tezduyar,  T.,  2003,  “Computation  of  Moving  Boundaries  and  Interfaces 
and  Stabilization  Parameters,”  Int.  J.  Numer.  Methods  Fluids,  43,  pp. 
555-575. 

[39]  John,  V.,  and  Knobloch,  P.,  2007,  “On  Spurious  Oscillations  at  Layers  Dimin¬ 
ishing  (Sold)  Methods  for  Convection-Diffusion  Equations:  Part  I — A  review,” 
Comput.  Methods  Appl.  Mech.  Eng.,  196(17-20),  pp.  2197-2215. 

[40]  Stein,  K.,  Tezduyar,  T.,  and  Benney,  R.,  2003,  “Mesh  Moving  Techniques  for 
Fluid-Structure  Interactions  With  Large  Displacements,”  J.  App.  Mech.,  70,  pp. 
58-63. 

[41]  Chung,  J.,  and  Hulbert,  G.  M.,  1993,  “A  Time  Integration  Algorithm  for  Struc¬ 
tural  Dynamics  With  Improved  Numerical  Dissipation:  The  Generalized-a 
Method,”  J.  App.  Mech.,  60,  pp.  371-375. 

[42]  Jansen,  K.  E.,  Whiting,  C.  H.,  and  Hulbert,  G.  M.,  1999,  “A  Generalized-a 
Method  for  Integrating  the  Filtered  Navier-Stokes  Equations  With  a  Stabilized 
Finite  Element  Method,”  Comput.  Methods  Appl.  Mech.  Eng.,  190,  pp. 
305-319. 

[43]  Saad,  Y.,  and  Schultz,  M.,  1986,  “GMRES:  A  Generalized  Minimal  Residual 
Algorithm  for  Solving  Non-Symmetric  Linear  Systems,”  SIAM  J.  Sci.  Comput., 
7,  pp.  856-869. 

[44]  Hughes,  T.  J.  R.,  and  Winget,  J.,  1980,  “Finite  Rotation  Effects  in  Numerical 
Integration  of  Rate  Constitutive  Equations  Arising  in  Large-Deformation  Ana¬ 
lysis,”  Int.  J.  Numer.  Methods  Eng.,  15,  pp.  1862-1867. 

[45]  Longo,  J.,  and  Stern,  F.,  2005,  “Uncertainty  Assessment  for  Towing  Tank  Tests 
With  Example  for  Surface  Combatant  DTMB  Model  5415,”  J.  Ship  Res.,  49, 
pp.  55-68. 

[46]  Garcia,  J.,  and  Onate,  E.,  2003,  “An  Unstructured  Finite  Element  Solver  for 
Ship  Hydrodynamics  Problems,”  J.  App.  Mech.,  70,  pp.  18-26. 

[47]  Longo,  J.,  Shao,  J.,  Irvine,  M.,  and  Stem,  F.,  2007,  “Phase-Averaged  PIV  for 
the  Nominal  Wake  of  a  Surface  Ship  in  Regular  Head  Waves,”  J.  Fluids  Eng., 
129,  pp.  524-541. 

[48]  McCormick,  M.,  2010,  Ocean  Engineering  Mechanics.  With  Applications, 
Cambridge  University  Press,  London. 

[49]  Bazilevs,  Y.,  and  Akkerman,  I.,  2010,  “Large  Eddy  Simulation  of  Turbulent 
Taylor-Couette  Flow  Using  Isogeometric  Analysis  and  the  Residual-Based 
Variational  Multiscale  Method,”  J.Comput.  Phys.,  229(9),  pp.  3402-3414. 


Journal  of  Applied  Mechanics 


JANUARY  201 2,  Vol.  79  /  010905-11 


